#This to obtain:
#Figure_A11a.jpg
#Figure_A11b.jpg
#Figure_A11c.jpg


library("ggplot2")
setwd("/Users/bgpopescu/")
#setwd("C:/Users/bogdanp/")


#Reading Polygons
HRV_adm0 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="HRV_adm0_wgs")

krajna <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                  layer="krajna_line_GCS")

HRV_adm3 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="HRV_adm3")
HRV_adm3<-st_simplify(HRV_adm3, dTolerance = 0.005)
HRV_adm3<-subset(HRV_adm3, select = c(ID_2))
HRV_adm3$ID_2

HRV_adm3_data = read_excel("./Dropbox/Legacies_Project/Analysis/data/merge.xlsx", sheet=1, col_names = TRUE, skip = 0)
HRV_adm3_data$ID_2


final_sp<-left_join(HRV_adm3, HRV_adm3_data, by = c("ID_2"="ID_2"))
final_sp$bosnia_border_mun
data<-subset(final_sp, bosnia_border_mun==0)
data$dist1 <- as.numeric( with (data,ifelse(data$treat==1, 1, -1)))
data$dist2 <-data$dist1*data$krajna6_distance




###########################
#PLACEBO: pct_no_water#
##########################

#Subset data 
placebo <-data[which(data$dist2>-100 & data$dist2<100),]


#Define range of locations for placebo border
sequence<-seq(from=-60, to=60, by=1)

#Create empty variables
p.value=rep(NA, length(sequence))
estimate=rep(NA, length(sequence))
confint1=rep(NA, length(sequence))
confint2=rep(NA, length(sequence))
#Loop over potential placebo locations of border
for(i in 1:length(sequence)){
  placebo$dist3=placebo$dist2+sequence[i]
  placebo$treat2=ifelse(placebo$dist3 > 0, placebo$treat <- 1, placebo$treat <- 0)
  model1<-lm(pct_no_water ~ treat2+
               POINT_X+POINT_Y+dist3+
               bfe1+bfe2+bfe3+bfe4+
               bfe5+bfe6+bfe7+
               zagreb_distance, data = placebo)
  #model1<-lm(pct_no_water~treat2+dist3+treat2*dist3, data=placebo)
  confint1[i]=confint(model1)[2,1]
  confint2[i]=confint(model1)[2,2]
  estimate[i]=coef(model1)[2]
}

placebo_data <- data.frame(sequence, estimate, confint1, confint2)

annotations <- data.frame(
  xpos = c(-Inf,-Inf,Inf,Inf),
  ypos =  c(-Inf, Inf,-Inf,Inf),
  annotateText = c("","← Closer to South \n  Bosnia",
                   "","Closer to North → \n  Slovenia and Hungary"),
  hjustvar = c(0,0,1,1) ,
  vjustvar = c(0,1,0,1)) #<- adjust


pct_no_water_placebo<-ggplot(data=placebo_data, aes(sequence, estimate))+
  geom_point(size=1.3)+
  
  geom_point(data = placebo_data[placebo_data$confint1 > 0.0 & placebo_data$sequence > -40,], color= "red", size = 2)+ 
  geom_errorbar(aes(ymin=confint1, ymax=confint2), size=0.2, width = 0)+
  geom_errorbar(data = placebo_data[placebo_data$confint1 > 0.0& placebo_data$sequence > -40,],
                aes(ymin=confint1, ymax=confint2), size=0.2, width = 0, color= "red")+
  geom_hline(yintercept = 0)+
  geom_vline(xintercept = 0)+
  theme_bw()+
  scale_x_continuous(name="Placebo border distance from actual border (km)", breaks = seq(-70, 70, 10)) +
  scale_y_continuous(name="Estimated effect")+
  ggtitle("Pct. Dwellings No Water, 2011") +
  theme(axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title=element_text(size=14),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))+
  geom_text(data=annotations,aes(x=xpos,y=ypos,hjust=hjustvar,vjust=vjustvar,label=annotateText))+
  geom_text(x=18,y=Inf, hjust=1,vjust=1,label="Actual Border", color="Red")


ggsave(pct_no_water_placebo, 
       file="./Dropbox/Legacies_Project/Paper/figures/Figure_A11a.jpg", 
       height=10, width=20, 
       units = "cm", dpi=300)


###########################
#PLACEBO: pov_rate_income#
##########################

#Subset data 
placebo <-data[which(data$dist2>-100 & data$dist2<100),]

#Define range of locations for placebo border
sequence=seq(from=-60, to=60, by=1)

#Create empty variables
p.value=rep(NA, length(sequence))
estimate=rep(NA, length(sequence))
confint1=rep(NA, length(sequence))
confint2=rep(NA, length(sequence))
#Loop over potential placebo locations of border
for(i in 1:length(sequence)){
  placebo$dist3=placebo$dist2+sequence[i]
  placebo$treat2=ifelse(placebo$dist3 > 0, placebo$treat <- 1, placebo$treat <- 0)
  model1<-lm(pov_rate_income ~ treat2+
               POINT_X+POINT_Y+dist3+
               bfe1+bfe2+bfe3+bfe4+
               bfe5+bfe6+bfe7+
               zagreb_distance, data = placebo)
  #model1<-lm(pov_rate_income~treat2+dist3+treat2*dist3, data=placebo)
  confint1[i]=confint(model1)[2,1]
  confint2[i]=confint(model1)[2,2]
  estimate[i]=coef(model1)[2]
}

placebo_data <- data.frame(sequence, estimate, confint1, confint2)

pov_rate_income_placebo<-ggplot(data=placebo_data, aes(sequence, estimate))+
  geom_point(size=1.3)+
  geom_point(data = placebo_data[placebo_data$confint1>0.0 & placebo_data$sequence<15 & placebo_data$sequence>-15,], color= "red", size = 2)+ 
  geom_errorbar(aes(ymin=confint1, ymax=confint2), size=0.2, width = 0)+
  geom_errorbar(data = placebo_data[placebo_data$confint1 > 0.0 & placebo_data$sequence<15 & placebo_data$sequence>-15,],
                aes(ymin=confint1, ymax=confint2), size=0.2, width = 0, color= "red")+
  geom_hline(yintercept = 0)+
  geom_vline(xintercept = 0)+
  theme_bw()+
  scale_x_continuous(name="Placebo border distance from actual border (km)", breaks = seq(-70, 70, 10)) +
  scale_y_continuous(name="Estimated effect")+
  ggtitle("Pct. People at Risk of Poverty, 2011") +
  theme(axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title=element_text(size=14),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))+
  geom_text(data=annotations,aes(x=xpos,y=ypos,hjust=hjustvar,vjust=vjustvar,label=annotateText))+
  geom_text(x=18,y=Inf, hjust=1,vjust=1,label="Actual Border", color="Red")

pov_rate_income_placebo
ggsave(pov_rate_income_placebo, 
       file="./Dropbox/Legacies_Project/Paper/figures/Figure_A11b.jpg", 
       height=10, width=20, 
       units = "cm", dpi=300)


####################
#PLACEBO: EDUCATION#
####################
#Subset data 
placebo <-data[which(data$dist2>-100 & data$dist2<100),]

#Define range of locations for placebo border
sequence=seq(from=-60, to=60, by=1)

#Create empty variables
p.value=rep(NA, length(sequence))
estimate=rep(NA, length(sequence))
confint1=rep(NA, length(sequence))
confint2=rep(NA, length(sequence))

#Loop over potential placebo locations of border
for(i in 1:length(sequence)){
  placebo$dist3=placebo$dist2+sequence[i]
  placebo$treat2=ifelse(placebo$dist3 > 0, placebo$treat <- 1, placebo$treat <- 0)
  model1<-lm(share_lesseduc_25 ~ treat2+
               POINT_X+POINT_Y+dist3+
               bfe1+bfe2+bfe3+bfe4+
               bfe5+bfe6+bfe7+
               zagreb_distance, data = placebo)
  confint1[i]=confint(model1)[2,1]
  confint2[i]=confint(model1)[2,2]
  estimate[i]=coef(model1)[2]
}

placebo_data <- data.frame(sequence, estimate, confint1, confint2)

edu_placebo<-ggplot(data=placebo_data, aes(sequence, estimate))+
  geom_point(size=1.3)+
  geom_point(data = placebo_data[placebo_data$confint1 > 0.0,], color= "red", size = 2)+ 
  geom_errorbar(aes(ymin=confint1, ymax=confint2), size=0.2, width = 0)+
  geom_errorbar(data = placebo_data[placebo_data$confint1 > 0.0,],
                aes(ymin=confint1, ymax=confint2), size=0.2, width = 0, color= "red")+
  geom_hline(yintercept = 0)+
  geom_vline(xintercept = 0)+
  theme_bw()+
  scale_x_continuous(name="Placebo border distance from actual border (km)", breaks = seq(-70, 70, 10)) +
  scale_y_continuous(name="Estimated effect")+
  ggtitle("Pct. Less Educated People, 2011") +
  theme(axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title=element_text(size=14),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))+
  geom_text(data=annotations,aes(x=xpos,y=ypos,hjust=hjustvar,vjust=vjustvar,label=annotateText))+
  geom_text(x=18,y=Inf, hjust=1,vjust=1,label="Actual Border", color="Red")

edu_placebo
ggsave(edu_placebo, 
       file="./Dropbox/Legacies_Project/Paper/figures/Figure_A11c.jpg", 
       height=10, width=20, 
       units = "cm", dpi=300)
